# categorical variables df
cat <- TrainTest %>% select(where(is.character))
# get the cardinality of a category
levels <- data.frame("vars" = names(cat),
                     "levels" = apply(cat, 2,
                     function(x) length(unique(x))) )
# detele rownames
rownames(levels) <- NULL

# sort by cardinality
levels_sorted <- levels %>% arrange(levels)

levels_sorted %>% kbl() %>% kable_minimal()

# select only categories with cardinality > 1 and < 2200
catvar <- levels %>%
  filter(levels < 3000, levels > 1) %>%
  select(vars)
# new dataframe
cat <- subset(cat, select = catvar$vars)

# select only categories with cardinality > 1 and < 2200
dum <- levels %>%
  filter(levels < 20, levels > 1) %>%
  select(vars)
# new dataframe
dummies <- subset(cat, select = dum$vars)

# numerical variables
num <- TrainTest %>% select(where(is.numeric))
# logical variables
logi <- TrainTest %>% select(where(is.logical))

# join three dataframes
traintest <- cbind(num, logi,cat) %>% as.data.table()
traintest <- traintest %>% mutate(across(where(is.logical), ~ as.numeric(.x)))

table(traintest$payment_type)

table(traintest$payment)

# quantity_group: mismas categorias que con la variable quantity
table(traintest$quantity_group)

# quantity_group: mismas categorias que con la variable quantity
table(traintest$quantity)

print_duplicates <- function(data) {
  duplicated_rows <- duplicated(data)
  duplicated_data <- data[duplicated_rows, ]
  duplicated_data %<>% as.data.frame()
  print(dim(duplicated_data)[1])
}

# duplicated rows
print_duplicates(traintest)

traintest <- traintest %>%
             select(-quantity_group) %>%
             select(-payment_type) %>%
             select(-installer)

# Function to impute zeros using missRanger method
impute <- function(data, columns, k, num_trees, impute = TRUE, seed = 0) {
  for (col in columns) {
    new_col <- paste0("NA_", col)
    # create new variable that says if the value was NA or not
    # if NA` -`> 1 ; else -> 0
    data[[new_col]] <- ifelse(data[[col]] == 0, 1, 0)
    data[[col]] <- ifelse(data[[col]] == 0, NA, data[[col]])
  }

  # now we impute values
  if (impute == TRUE) {
    data <- data %>%
            # create variables that tell us if the row was NA or not
            mutate(NA_public_meeting = ifelse(is.na(public_meeting),0,1)) %>%
            mutate(NA_permit = ifelse(is.na(permit),0,1))
    data_imp <- missRanger(data,
                           pmm.k = 5,
                           seed = 0,
                           maxiter = 100)
    #data_imp <- kNN(data)#,k = floor(sqrt(dim(data)[1])))

  }
  return(data_imp)
}

# Function to impute zeros using KNN method
impute <- function(data, columns, k, num_trees, impute = TRUE, seed = 0) {
  # now we impute values
  if (impute == TRUE) {
    data <- data %>%
            # create variables that tell us if the row was NA or not
            mutate(NA_public_meeting = ifelse(is.na(public_meeting),0,1)) %>%
            mutate(NA_permit = ifelse(is.na(permit),0,1))
    data_imp <- kNN(data)#,k = floor(sqrt(dim(data)[1])))

  }
  return(data_imp)
}

# Function to impute zeros using MICE method
impute <- function(data, columns, k, num_trees, impute = TRUE, seed = 0) {
  # now we impute values
  if (impute == TRUE) {
    data <- data %>%
            # create variables that tell us if the row was NA or not
            mutate(NA_public_meeting = ifelse(is.na(public_meeting),0,1)) %>%
            mutate(NA_permit = ifelse(is.na(permit),0,1))
    data_imp <- mice(data, m = 20, method = "pmm", seed = 0)
    data_imp <- complete(data_imp)
  }
  return(data_imp)
}

# Funcion para la imputacion de ceros con missRanger
impute <- function(data, columns, k, num_trees, impute = TRUE, seed = 0) {
  # now we impute values
  if (impute == TRUE) {
    data <- data %>%
            # create variables that tell us if the row was NA or not
            mutate(NA_public_meeting = ifelse(is.na(public_meeting),0,1)) %>%
            mutate(NA_permit = ifelse(is.na(permit),0,1))
    data_imp <- missRanger(data,
                           pmm.k = 5,
                           seed = 0,
                           maxiter = 100)
  }
  return(data_imp)
}

# impute data
traintestimp <- impute(traintest,
                             columns = c('construction_year', 'district_code', 'longitude','gps_height'),
                             k = 5,
                             num_trees = 100)

# cleanning some garbage present in ward and funder to reduce categories
clean_text <- function(text) {
  stri_trans_tolower(
    stri_replace_all_regex(
      text,
      pattern = "[ +\\p{Punct}()]",
      replacement = "")
    )
}

traintestimp <- traintestimp %>%
        #population por bomb number * ward
        #mutate(fe_popnum_ward = (sum(population) / n.()), .by = ward ) %>%
        mutate(ward = clean_text(ward)) %>%
        mutate(funder = clean_text(funder)) %>%
        mutate(scheme_name = clean_text(scheme_name)) %>%
        mutate(ward = as.character(fct_lump(ward, n = 200))) %>%
        mutate(funder = as.character(fct_lump(funder, n = 200))) %>%
        mutate(scheme_name = as.character(fct_lump(scheme_name, n = 200))) %>%
        #geographic distance
        mutate(fe_dist = distGeo(as.matrix(cbind(longitude, latitude)), c(0, 0))) %>%
        # amount of water per population
        mutate(fe_cant_agua = ifelse(population == 0, 0, round(amount_tsh / population, 3))) %>%
        # age of the water pump
        mutate(fe_age = 2014 - construction_year) %>%
        mutate(fe_days = as.numeric(as.Date("2014-01-01") - as.Date(TrainTest$date_recorded))) %>%
        # get the year and the month
        mutate(fe_dryear = year(TrainTest$date_recorded)) %>%
        mutate(fe_drmonth = month(TrainTest$date_recorded)) %>%
        # difference between recorded date and construction date
        mutate(fe_drdiff = fe_dryear-construction_year) %>%
        as.data.table()

cat_cols <- names(traintestimp %>% select(where(is.character)))
cat_cols

for (cat_col in cat_cols) {
  # calculate for every value and then convert it to numeric
  traintestimp[, paste0("fe_", cat_col) := as.numeric(.N), by = cat_col] }
# delete old cat cols
for (cat_col in cat_cols) {
  traintestimp[, paste(cat_col) := NULL] }

datos_dummy <- dummy_cols(dummies) %>% select(where(is.numeric))
traintestimp <- cbind(traintestimp, datos_dummy)

numeric_vars <- names(traintestimp %>%  select(where(is.numeric)))
length(numeric_vars)

# Función para estandarizar las variables
standardize <- function(x) {
  (x - mean(x)) / sd(x) }

# Escalar las variables numéricas
traintestimp[, (numeric_vars) := lapply(.SD, standardize), .SDcols = numeric_vars]

traintestimp %<>% select(-num_private) %>%
                  select(-NA_permit) %>%
                  select(-NA_public_meeting)  %>%
                  select(-fe_dryear)

# create idx column to separate
traintestimp %<>%
  mutate( idx = 1:nrow(traintestimp)) %>%
  as.data.table()
# separate train and test using idx.
xtrain <- traintestimp %>% as.data.table() %>%
  filter( idx <= nrow(train)) %>%
  select(-idx) %>%
  #mutate( status_group = trainlab$status_group) %>% # we add the target
  as.data.table()

xtest <- traintestimp %>% as.data.table() %>%
  filter( idx > nrow(train)) %>%
  select(-idx) %>%
  as.data.table()

vector_status_group <- trainlab$status_group
vector_status_group <- ifelse(vector_status_group == "functional", 0,
                              ifelse(vector_status_group == "functional needs repair", 1,
                              2))

vector_status_group <- as.factor(vector_status_group)
x_train <- xtrain %>% mutate(vector_status_group=vector_status_group) %>%   as.data.frame()
x_test <- as.data.frame(xtest)

formula <- as.formula(vector_status_group ~ .)

table(vector_status_group)

# completely unbalanced in class 1

# Apply SMOTE
set.seed(0)

data <- x_train %>% mutate(vector_status_group = as.factor(vector_status_group))

# to increase the number of samples of the minoritary class
smote_data <- SMOTE(vector_status_group ~ ., data)#, perc.over = 100,perc.under = 500)

# Verify if smote worked correctly
table(smote_data$vector_status_group)

# Definir los parámetros
num_trees <- 100  # Número de árboles en el bosque
mtry <- 3  # Número de variables predictoras seleccionadas en cada división
n_folds <- 5  # Número de folds en la validación cruzada

# Crear un contenedor para almacenar los resultados de cada fold
accuracy_results <- numeric(n_folds)

# Realizar la validación cruzada
set.seed(42)  # Fijar la semilla para reproducibilidad

#x_train <- xtrain %>% mutate(vector_status_group=vector_status_group) %>%   as.data.frame()
# Obtener los índices de los folds
folds <- cut(seq(1, nrow(x_train)), breaks = n_folds, labels = FALSE)


rf_accuracy <- numeric(n_folds)
logistic_accuracy <- numeric(n_folds)
xgb_accuracy <- numeric(n_folds)
lgb_accuracy <- numeric(n_folds)
svm_accuracy <- numeric(n_folds)

# Iterar sobre cada fold
for (i in 1:n_folds) {
  # Dividir los datos en conjuntos de entrenamiento y prueba según el fold actual
  train_data <- x_train[folds != i, ]
  test_data <- x_train[folds == i, ]

  vector_train <- train_data %>% select(vector_status_group) %>% pull()
  train_matrix <- train_data %>% select(-vector_status_group) %>% as.matrix()

  vector_test <- test_data %>% select(vector_status_group) %>% pull()
  test_matrix <- test_data %>% select(-vector_status_group) %>% as.matrix()


  params <- list(
  objective = "multiclass",
  num_class = 3,
  boosting_type = "gbdt",
  metric = "multi_logloss")
  # learning_rate = 0.02,
  # max_depth = 15,
  # num_threads = 4)
  # Convertir los datos en formato xgb.DMatrix
  lgb_train <- lgb.Dataset(data = train_matrix, label = as.numeric(vector_train) - 1)
  lgb_test <- lgb.Dataset(data = test_matrix, label = as.numeric(vector_test) - 1)
  # Ajustar el modelo de lightgbm en el conjunto de entrenamiento
  capture.output({ lgbm <- lgb.train(data = lgb_train, nrounds = 600, valids = list(test = lgb_test), early_stopping_rounds = 10, params = params) }, file = NULL) #params = params,)
  pred_lgb <- predict(lgbm, test_matrix)
  pred_lgb_class <- max.col(matrix(pred_lgb, ncol = 3, byrow = TRUE)) - 1
  lgb_accuracy[i] <- sum(pred_lgb_class == as.numeric(vector_test) - 1) / length(pred_lgb_class)

  # Ajustar el modelo de Random Forest en el conjunto de entrenamiento
  rf <- ranger(formula, importance = 'impurity', data = train_data, num.trees = num_trees, mtry = mtry, verbose = FALSE, seed = 42, probability = TRUE)
  predrf <- predict(rf, data = test_data)$predictions
  rf_accuracy[i] <- sum(max.col(predrf) == as.numeric(test_data$vector_status_group)) / length(test_data$vector_status_group)

  # Ajustar el modelo de regresión logística multinomial en el conjunto de entrenamiento
  glm <- suppressMessages(multinom(formula, data = train_data, MaxNWts = 10000, MaxIter = 100))
  predglm <- predict(glm, newdata = test_data, type = "class")
  logistic_accuracy[i] <- sum(predglm == test_data$vector_status_group) / length(test_data$vector_status_group)

  # Ajustar el modelo de xgboost en el conjunto de entrenamiento
  params <- list(
    objective = "multi:softmax",
    num_class = 3,
    eval_metric = "mlogloss")

  # Entrenar el modelo XGBoost
  xgb_train <- xgb.DMatrix(data = train_matrix, label = as.numeric(vector_train) - 1)
  xgb_test <- xgb.DMatrix(data = test_matrix, label = as.numeric(vector_test) - 1)

  xgb_model <- xgb.train(params = params, data = xgb_train, nrounds = 600, nthread = 4)
  predxgb <- predict(xgb_model, xgb_test)
  xgb_accuracy[i] <- sum(predxgb == as.numeric(vector_test) - 1) / length(predxgb)

  # Ajustar el modelo de svm en el conjunto de entrenamiento
  train_data_svm <- data.frame(train_matrix, y = as.factor(as.numeric(vector_train) - 1))
  test_data_svm <- data.frame(test_matrix, y = as.factor(as.numeric(vector_test) - 1))
  svm_model <- svm(y ~ ., data = train_data_svm, method = "svmRadial", kernel = "radial", cost = 10, scale = TRUE)
  pred_svm <- predict(svm_model, test_data_svm)
  svm_accuracy[i] <- sum(pred_svm == as.numeric(vector_test) - 1) / length(pred_svm)
}

# Define the parameter grid to explore
grid <- expand.grid(
  mtry = c(6,7,8),                 # Values for the number of predictor variables considered at each split
  min.node.size = c(4, 5, 6),        # Values for the minimum node size
  splitrule = c("gini") # Values for the splitting rule
)

# Define the control function for hyperparameter search
ctrl <- trainControl(
  method = "cv",    # Use cross-validation
  number = 5,       # Number of cross-validation folds
  search = "grid"   # Perform grid search on the defined parameter grid
)

# Perform the hyperparameter search and train the model
rf_model <- train(
  formula,
  data = x_train,
  method = "ranger",
  trControl = ctrl,
  tuneGrid = grid
)

# Show the best selected parameters
print(rf_model$bestTune)

train_matrix <- x_train %>% select(-vector_status_group) %>%  as.matrix()
test_matrix <- xtest %>% as.matrix()

# Convert the data to DMatrix format
dtrain <- xgb.DMatrix(data = train_matrix, label = as.numeric(vector_status_group) - 1)
dtest <- xgb.DMatrix(data = test_matrix)

pb <- progress_bar$new(total = nrow(params_grid), format = "[:bar] :percent ETA: :eta")

# Define the hyperparameters for the grid search
params_grid <- expand.grid(
  max_depth = c(6, 8, 10),
  eta = c(0.01, 0.1),
  gamma = c(0, 1),
  colsample_bytree = c(0.5, 0.75, 1),
  min_child_weight = c(1, 5),
  subsample = c(0.5, 0.8),
  objective = "multi:softprob",
  num_class = length(unique(vector_status_group)) )


best_acc <- 0
best_params <- list()

pb <- progress_bar$new(total = nrow(params_grid), format = "[:bar] :percent ETA: :eta")

for (i in 1:nrow(params_grid)) {
  params <- list(
    max_depth = params_grid$max_depth[i],
    eta = params_grid$eta[i],
    gamma = params_grid$gamma[i],
    colsample_bytree = params_grid$colsample_bytree[i],
    min_child_weight = params_grid$min_child_weight[i],
    subsample = params_grid$subsample[i],
    objective = "multi:softprob",
    num_class = length(unique(vector_status_group))
  )
  # Perform grid search for hyperparameter tuning
  xgb_model <- xgb.cv(
    params = params,
    data = dtrain,
    nfold = 5,
    nrounds = 200,
    early_stopping_rounds = 20,
    verbose = 0
  )

   if (1 - tail(xgb_model$evaluation_log$test_mlogloss_mean, 1) > best_acc) {
    best_acc <- 1 - tail(xgb_model$evaluation_log$test_mlogloss_mean, 1)
    best_params <- params
  }
  # progress bar to show the %
  pb$tick()
}

print(best_params)

split_index <- createDataPartition(x_train$vector_status_group, p = 0.7, list = FALSE)
x_train_base <- x_train[split_index, ]
x_train_meta <- x_train[-split_index, ]

# Convertir x_train_base y x_train_meta a matrices para XGBoost
y_train_base <- vector_status_group[split_index]
y_train_meta <- vector_status_group[-split_index]

x_train_base_matrix <- x_train_base %>% select(-vector_status_group) %>%  as.matrix()
x_train_meta_matrix <- x_train_meta %>% select(-vector_status_group) %>%  as.matrix()

# Ranger model
ranger_model <- ranger(formula,
                       data = x_train_base, num.trees = 500,
                       mtry = 7, importance = 'impurity',
                       min.node.size = 5, splitrule = 'gini',
                       verbose = FALSE, seed = 42, probability = FALSE)

# Predecir con el modelo Ranger
ranger_predictions_meta <- predict(ranger_model, x_train_meta, type = "response")$predictions

# Xgboost model
params <- list(max_depth = 15, eta = 0.02,
               gamma = 0, colsample_bytree  = 0.3,
               colsample_bylevel = 0.9, colsample_bynode  = 0.7,
               min_child_weight = 1, subsample = 0.8,
               objective = "multi:softmax", num_class = 3)

# Convert the data to DMatrix format
dtrain_base <- xgb.DMatrix(data = x_train_base_matrix, label = as.numeric(y_train_base) - 1)
# train model
xgb_model <- xgb.train(params = params, data = dtrain_base, nrounds = 600)

# Predecir con el modelo XGBoost
xgb_predictions_meta <- predict(xgb_model, x_train_meta_matrix)

# Crear un data.frame con las predicciones de ambos modelos
combined_predictions_meta <- data.frame(ranger = ranger_predictions_meta, xgboost = xgb_predictions_meta)

# Entrenar el modelo meta
#meta_model <- randomForest(y ~ ., data = cbind(y = x_train_meta$vector_status_group, combined_predictions_meta), ntree = 500)
library(nnet)
meta_model <- multinom(y ~ ., data = cbind(y = x_train_meta$vector_status_group, combined_predictions_meta))

# Predecir con el modelo Ranger
ranger_predictions_test <- predict(ranger_model, x_test, type = "response")$predictions

# Predecir con el modelo XGBoost
test_matrix <- xtest %>% as.matrix()
xgb_predictions_test <- predict(xgb_model, test_matrix)

# Crear un data.frame con las predicciones de ambos modelos
combined_predictions_test <- data.frame(ranger = ranger_predictions_test, xgboost = xgb_predictions_test)
# Stacking de los modelos Ranger y XGBoost
stacked_predictions <- predict(meta_model, combined_predictions_test)

# Make predictions on the test dataset
stacked_predictions <- ifelse(stacked_predictions == 0, "functional", ifelse(stacked_predictions == 1, "functional needs repair", "non functional"))

# Submit predictions
stacked_predictions <- data.table(id = test$id, status_group = stacked_predictions)
fwrite(stacked_predictions, file = "./submissions/stacking")